% This function generate the sample sequences from the signal and the
% specified MWC system

function [SignalSampleSequences, Sample_t_axis] = MWC_SubNyquistSample(Signal,MWC, SignalAmpAnalog, noise_unscaled, Analog_TimeAxis,SNR)

% The procedure of generating the lowrate samples, in numerical simulation,
% consists of the following stages:
% 1. generate the input signal at the Nyquist rate
% 2. generate the sequences z[n], by appropriate mixing and filtering.
% These sequences are decimated to rate f_s.
% 3. Let q=floor(Tp/Ts); A=S*F*D ; if q>1 then update A according to the
% expansion factor q
% 4. Mix the sequences z[n] according to the matrix A, and decimate the
% results to 1/Ts rate

% 1: generate the signal at the Nyquist rate (or higher)

% We would like to avoid noninteger decimations, therefore we generate the
% signal on a time-axis that is at least at the Nyquist resolution,
% and that allows for integer decimations in the sequel

%factor to resemble the signal to analog one 
AnalogFactorRes = 5;



fmax = getfield(Signal.Structure, 'fmax');
fnyq = 2*fmax;
M = getfield(MWC, 'M');
m = getfield(MWC, 'm');
fp = getfield(MWC, 'fp');
Ts = getfield(MWC, 'Ts');
fs = getfield(MWC, 'fs');
q= getfield(MWC, 'q');
decfactor = ceil(fnyq*AnalogFactorRes/fs); % this is the decimation ratio from the Nyquist-rate to fp*q
fgrid = decfactor*fs;  % fgrid>=fnyq, and both fgrid/fp and fgrid/fs are integers
% decfactor = ceil(fnyq/fp); % this is the decimation ratio from the Nyquist-rate to fp
% fgrid = decfactor*fp;  % fgrid>=fnyq, and both fgrid/fp and fgrid/fs are integers


% % get the required number of samples
% NumSamples = 450*q ; % 999 Understand the mining of this number   %%str2double(get(handles.editNumSamples,'String'));
% Npts = NumSamples*decfactor*2;  % we take a factor of 2.. just to be safe
% 
% % generate signal at analog like density
% [SignalAmpAnalog, NoiseAmpAnalog, Analog_TimeAxis] = GenerateSignalValues(Signal,1/fgrid,Npts);

%% adapt noise to desired SNR level
NoiseEnergy = norm(noise_unscaled)^2;
SignalEnergy = norm(SignalAmpAnalog)^2;
CurrentSNR = SignalEnergy/NoiseEnergy;
SNR_dB = SNR;
SNR_val = 10^(SNR_dB/10);
NoiseAmpAnalog = noise_unscaled*sqrt(CurrentSNR/SNR_val);  %%999 considre moving this stage to later phase and allow scan through different SNRs
%%
% figure(1)
% plot(SignalAmpAnalog);
% hold all;
% plot(NoiseAmpAnalog);
% hold off;
% legend('signal','noise');
% title(['SNR of ', num2str(SNR)])
% xlabel('time[t]');
% ylabel('Amplitude')

%%
Sig = SignalAmpAnalog; %ETGAR
% Sig = SignalAmpAnalog + NoiseAmpAnalog;

%% Plot Signal and Noise in the frequency domain
% Tanalog = Analog_TimeAxis(2)-Analog_TimeAxis(1);% Sample time
% title_1='Double-Sided Amplitude Spectrum of the noisy signal(t) - pre to any processing';
% legendText={'Signal';'Noise'}
% % display_at_freq_domain( 1/Tanalog,Sig,fmax,title_1,legendText);
% display_at_freq_domain( 1/Tanalog,[SignalAmpAnalog;NoiseAmpAnalog],fmax,title_1,legendText,'b');
%%

% 2: Generate z[n]

L = MWC.L;

if (mod(L,2)==0)
    L=L+1;
end   % L is the number of specturm slices, namely the length of z[n]

% prepare a polyphase filter H(f)   
% Analog filter (polyphase for speed)
% dec = ceil( fgrid/(fs-fp));
dec = fgrid/fs;
h_analog = fir1(5000, 1/dec);
lenh = length(h_analog);
polylen = ceil(lenh / dec);
h_pad = [h_analog zeros(1,polylen*dec-lenh)];
polyh = reshape(h_pad,dec,polylen);

%% DEBUG - ANALYZE ANLALOG FILTER PROPERTIES
% fvtool(h_analog,1);
%%
%% 3: Analog Mixing
SignPatterns = getfield(MWC, 'SignPatterns');

MixedSigSequences = zeros(m,length(Analog_TimeAxis));
for channel=1:m
    MixedSigSequences(channel,:) = MixSignal(Sig,Analog_TimeAxis,SignPatterns(channel,:),1/fp);
end

%% Display the mixed signal at the frequency domain
% title_2='Double-Sided Amplitude Spectrum of the noisy signal(t) -after the mixing';
% display_at_freq_domain(1/Tanalog,MixedSigSequences(4,:),fmax,title_2);

%% Analog low-pass filtering and actual sampling

SignalSampleSequences = [];
for channel = 1:m
    SignalSequence =  MixedSigSequences(channel,:);

    [tempSampledSignal, Sample_t_axis] = PolyFilterDecimate(SignalSequence,Analog_TimeAxis,polyh,lenh);
    SampleLength = length(Sample_t_axis);

    SignalSampleSequences(channel, 1:SampleLength) = tempSampledSignal;
    
end
%% Bypass Analog MWC
A  = computeSensingMatrix( MWC);
Sig = downsample(Sig,5);
size(A)
A = repmat(A,1,ceil(length(Sig)/255));
A = A(:,1:length(Sig));
size(Sig)
size(A)
Sig = repmat(Sig',1,length(Sig));
size(Sig)
SignalSampleSequences = A*Sig.';


%% DISPLAY SIGNAL AND NOISE AFTER FOLDING THE FREQUENCY DOMAIN INTO THE BASEBAND
% title_3='Double-Sided Amplitude Spectrum of the noisy signal(t) -after the folding';
% display_at_freq_domain( fs,SignalSampleSequences(4,:),fs/2,title_3);
end

